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Numerical solution of the time -dependent compressible 
Navier-Stokes equations in inlet regions 
By L. E. OLSON 

United States Army Air Mobility Research and Development Laboratory, 

Ames Research Center 
P. R. McGOWAN 

Computer Sciences Corporation, Field Services Division, Mountain View, 

California 94043 
AND R. W. MacCORMACK 

Ames Research Center, National Aeronautics and Space Administration, 
Moffett Field, California 94035 

This paper presents the results of a study to determine the effects of 
compressibility on the viscous flow through channels that have straight, 
parallel walls. Two channel configurations are considered, the flow between 
two semi-infinite flat plates with uniform flow prescribed at the inlet plane 
and a cascade of semi- infinite flat plates with uniform flow introduced 
upstream. Tlie flow field is modeled by using the time-dependent, compressible 
Navier-Stokes equations. Time -dependent solutions are obtained by using an 
explicit finite-difference technique which advances the pressure on near-field 
subsonic boundaries such that accurate steady-state solutions are obtained. 

The steady-state results at Reynolds numbers 20 and 150 are presented for Mach 
numbers between 0.09 and 0.36 and compared with the incompressible solutions 
of previous studies. 
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1 . Introduction 

The flow in the inlet region of a straight channel brings with it many 
of the features typical of viscous flow, and thus provides an interesting test 
case for developing new analytical and numerical methods to be used in the 
study of more practical flows. As a result, there is an extensive list of 
references which address this problem, beginning with the original work of 
Boussinesq (1891). Analytical studies by Schlichting (1960), Blankenship and 
Chung (1967), Van Dyke (1970), Wilson (1971), and Kapila, Ludford, and 
Olunloyo (1970) represent some of the more current publications where asymptotic 
solutions for large Reynolds number have been constructed. With the exception 
of the work by Blankenship and Cheng, all are restricted to incompressible 
flow. Additional investigators have approached the problem numerically. 

Among those using the incompressible Navier-Stokes equations are Gillis and 
Brandt (1964), Wang and Longwell (1964), Freidman, Gillis, and Liron (1968), 
Gosman (1969), and McDonald, Denny, Mills (1972), and Morihara and Cheng 
(1973) . Taylor and Ndefo (1970) solve the time-dependent incompressible 
equations of motion. These works, and others, form the basis for a fundamental 
understanding of incompressible inlet flows. To the authors' knowledge, how- 
ever, no solutions have been obtained for this type of problem using the 
unsteady, compressible Navier-Stokes equations. 

Three basic models are generally considered in the literature: uniform 

flow into a straight channel, irrotational flow into a straight channel, and 
a cascade of semi-infinite flat plates immersed in a uniform flow. The infinite 
cascade most accurately simulates a physically realizable flow, and is con- 
sidered in this report, as is the uniform entry flow which Van Dyke (1970) 
suggests may be attainable through the use of a variable’ porosity mesh. In 
flows such as these, it is important that particular attention be given to the 
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boundary conditions on near-field permeable boundaries; that is, boundaries 
where the flow either enters or leaves the computational region. This situa- 
tion is relieved to some degree for incompressible flow, where it is possible 
to eliminate the pressure from the governing equations by a transformation 
into the stream function -vorticity plane. The resulting equations, when 
coupled with the appropriate boundary conditions, can then be used to solve for 
the velocity field without reference to any boijndary condition on the pres- 
sure. The pressure field can then be computed by integrating the momentum 
equations with a priori knowledge of the pressure required at only one point 
in the flow field. 

For the compressible case, it is not possible to uncouple the pressure 
and, as a result, a boundary condition on this variable is also required. In 
addition, with unsteady calculations the pressure conditions become a function 
of time, which must be generated as part of the solution. In the present 
paper, a boundary condition which permits the pressure on an upstream subsonic 
inflow boundary to become a function of time is combined with the explicit 
finite-difference technique of MacCormack (1970) to solve the compressible 
time -dependent Navier-Stokes equations. This boundary condition does not 
require a priori knowledge of the static pressure on the inflow boundary, 
although the analysis does assume that boundary conditions on the two velocity 
components and density at this location are specified. Steady-state solutions 
with low subsonic Mach numbers are compared with the incompressible results of 
previous authors. Solutions at Mach numbers where compressibility effects 
become important are also presented. 

2. Differential equations and numerical method 

The Navier-Stokes equations, in two dimensions, for a time-dependent 
compressible flow may be written in vector form as 
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In the above equations, p represents the static pressure, p is the density, 
T is the static temperature, u is the component of velocity in the x- 
direction, and v is the component of velocity in the y-direction* 

The additional equations used to close the system are the equation of 


state. 
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p = pRT, 


and the total energy per unit volume. 


e = pe + p 




where c is the internal energy per unit mass, which is related to the static 


temperature through 


Sutherland's law. 
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was used for the viscosity, X - -(2/3)y, and the Prandtl number was taken to 
be 0.72. The ratio of specific heats, y, was set equal to 1.4. 

The split predictor-corrector difference method devised by MacCormack 
C1970) to solve equation (1) was used, and can be defined as follows: 
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This essentially consists of two applications of a modified Euler predictor- 
corrector pair. The first application, equations (2), ignores the x- 
derivatives in equation (1)] whereas, the second ignores the /-derivatives. 
This numerical computation uses data obtained from dqta in the t = t^^ plane, 
see figure 1, to advance the solution at x^, y^ by an amount At. 

If Ly^(At) is defined as an operator which represents the operations 
performed in equations (2) and L^(At) is the corresponding operator for equa- 
tions (3), then the solution of equation (1) can be advanced in time, with 
second-order accuracy, by the following sequence: 

= Ly(At)Lx(At)Lx(At)Ly(4t)U^ . (4) 


Equation (4) holds for Ay = Ax; however, significant advantages result 
when a modification of this technique is used in problems in which it is 
desirable to have Ay Ax. In this case, larger time steps and closer match- 
ing of the numerical and physical domains of dependence is achieved as follows: 
The solution of equation (1) is advanced by an amount At^^. with 


^ ^M/2Ly(Aty)Lx(Atx)iT^^^Ly(Aty)u",, 




where M is the smallest even integer greater than At„/At,y., and ir 


M/2 


denotes 


M/2 products of the Ly operator. For the unequal spatial step sizes, 
stability considerations require that 
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where u and c are the local velocity and speed of sound, respectively. 

A rectangular grid which grows geometrically in both the x- and the 
/-direction is shown in figure 1, along with the coordinate system. This 


<Fig' 
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arrangement is ideally suited to the inlet flow problem, which requires a fine 
mesh in the immediate vicinity of the leading edge but permits a more coarse 
mesh in the vicinity of the centerline, downstream, and upstream. For the 
geometric grid, the mesh spacing is given by 

AXi = (1 + n^)Ax^_^ = Cl + n^)^‘^Ax, 

and 

= (1 + 

In general, n was not equal to n and, in the case of the cascade, the r\^ 
X y 

used upstream of the leading edge was not necessarily equal to the used 

to the right of the leading edge. 

The Ax and Ay required for good spatial resolution is estimated by 

assuming that significant terms in the differential equations are of the same 

order of magnitude. Consider the transport term 9(pvu)/3y and the viscous 

stress term 3 [y (9u/9y)]/9y in the x-momenturn equation. Linerization and 

differencing gives, for n = 0 = 0, 
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If these terms are of the same order of magnitude, then one would expect the 
coefficients also to be of the same order of magnitude; that is, 
pv/2Ay ~ u/2CAy)^. This condition, when expressed in terms of mesh Reynolds 
number, implies 



pvAy ^ 2 
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or, upon introducing the Reynolds number. Re - 2hu^/v^, 
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where h is the half-channel width. The results of Gillis and Brandt (1964) 
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indicate that, for a Reynolds number of 200, jv/u^| « 1/6, which indicates 
Ay/h z 0.08, is sufficient for good spatial resolution at Reynolds numbers 
at or below 200. 

Equating linearized and differenced terms 3(pu^)/9x and 8(pvuV9u from 
the y-momentum equations yields Ax - (u/v)Ay or, again using the results of 
Gill is and Brandt, Ax 2 6Ay. In the present calculations, Ay/h was 0.036 or 
less and Ax was fixed at 2Ay, both of which are conservative according to 
the above estimates for good spatial resolution. 

With Ax determined, the number of mesh points to be used from the inlet 
plane to the downstream exit plane (N) specified, and the distance from the 
inlet plane to the exit plane Cx ) specified, the geometric scale factor to 
the right of the inlet plane (hx) is computed using the equation 


e 

h 


AX 

h 


(1 - - 1 


The scale factors to the left of the inlet and in the y-direction are computed 
using similar formulae. These equations center the smallest cell in the 
numerical grid (which has dimensions Ax, Ay) on the leading edge. 


3. Initial conditions, boundary conditions, and convergence criteria 

The initial conditions are essentially that of a uniform flow, which in 
terms of the present notation, is u = u^, p = p^, p = p^, and T = T^ for 
all X and y. The posing of boundary conditions in any numerical calculation 
requires particular attention. This is particularly true in subsonic flow. 

As noted by Moretti (1969), the boundary of a subsonic computational region with 
a permeable boundary on which flow variables are held constant is an ill- 
posed problem. Certainly this is the case for the time-dependent inlet flow. 
Perturbation waves created in the interior of the flow accumulate at this type 
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o£ boundary. Moretti observed this characteristic in Laval nozzle computa- 
tions, and similar behavior has been observed in the present calculations for 
the case of the uniform inlet when ‘all flow variables at the inlet plane were 
held fixed. Thus, it is necessary to permit the boundary conditions to 
become functions of time without altering the fundamental character of the 
problem being studied. 

First, consider the boundary conditions for the case of uniform entry. 

The boundary conditions at the inlet plane specify uniform flow and constant 

density, which in the present notation is u. . = u , v. . = 0, and p.. = p . 

•' ^ ij ^ 13 ij 

The numerical computation of Gill is and Brandt (1964) and the analytic work 
of Van Dyke (1970) indicate that a pressure gradient across the channel inlet 
plane is required to maintain the prescribed uniform velocity profiles. There- 
fore, the pressure distribution across the face of the inlet is permitted to 
adjust by linear extrapolation from the first two interior points to the inlet 
boundary after each time step. This boundary condition is then written as 

Axi 

n+1 n , n n . ^ 

^ij ~ ^2j ^3j Ax2 

The inlet temperature distribution is advanced by using the equation of state. 

For the boundary conditions at the plate, a fictitious row of cells 
below the plate were used, together with the reflection technique on u and v. 
The pressure in the fictitious cell was updated during each time step by 
using, as the compatibility condition, the y-raomentum equation evaluated at the 
surface 

9p _ u 9^u _y_ ^ dj_ 4^ p9^v 

9y 3 3x3y 2T 3y 3y 3 ^^2 

For an adiabatic wall, the coirplete set of boundary conditions then becomes 


“il = -“i2> ''il = -^2' hi = h2- 


and 
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Symmetry boundary conditions were applied on the channel centerline by using 
mirror cells on either side of the channel centerline. 

The downstream boundary values on u, v, and p were updated by linear 
extrapolation from the two upstream points. This is equivalent to forcing 
9^u/3x^ = 9^v/9x^ = 3^p/9x^ = 0. The pressure on this boundary, which was 
initially uniform, was held fixed at the initial value. The equation of state 
was again used to update the temperature. 

The distance to the parabolic regime as determined for incompressible inlet 
flow by Gillis and Brandt (1964) was used to determine the approximate location 
of the downstream boundary. The basic requirement is that the downstream 
boundary be sufficiently far removed from the inlet to permit the pressure 
distribution to become essentially uniform across the channel and, therefore, 
consistent with the specified downstream boundary condition of uniform pressure 
at the exit plane. In all of the cases presented in this paper, including 
both the uniform inlet and cascade calculations, the normalized static pressure 
distribution, p/p^u^, at the first interior plane upstream of the downstream 
boundary, was uniform to within 1% for the converged solutions. 

The boundary conditions for the cascade are the same as those used for 
the uniform inlet flow, except ahead of the plate where symmetry requires that 
V = 9u/9y = 9p/9y = 3p/3y = 0 on the lower boundary. The distance from the 
cascade inlet to the upstream boundary was located by specifying that the con- 
verged pressure distribution Cp/p^u^) over the upstream plane be uniform to 
within 0.5%. 

Equation (1), together with the appropriate boundary condition, was 
solved for the uniform inlet flow and the cascade at various combinations of 
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Reynolds number and Mach number, M^, ranging from 20 to 150 and 0.089 to 
0.364, respectively. For both the inlet and the cascade, convergence to the 
steady-state solution was not monotonic. Rather, exponentially damped sinu- 
soidal oscillations were observed. Disturbances created in the startup process 
are at least partially reflected at the channel inlet plane and the downstream 
exit plane and as a result persist in the interior of the flow field until 
damped by transmission at the boundaries, viscous dissipation, and heat- 
conduction effects. Although the somewhat arbitrary method used to update the 
pressure on the upstream boundary precludes precise comparison, the observed 
frequencies for both the uniform inlet and cascade flows are in rough agree- 
ment with the fundamental resonant frequencies of channels which have one end 
closed the other end open to the atmosphere (Malecki 1969) . In addition, 
the dependence of the rate of dissipation on frequency showed approximate 
agreement with the analysis of Malecki, which indicates that radiation, 
viscous dissipation, and dissipation due to heat condution from a closed-end 
channel are proportional to the square of the resonant frequency. The solu- 
tions with Re = 150 were assumed to be converged to the steady-state solution 
when An/o u^ 1 0.02, where Ap is defined as the amplitude of the pressure 
fluctuation on the flow centerline at the upstream boundary. The Re = 20 
solutions, which converged considerably faster than the Re = 150 cases, were 
assumed converged when Ap/p^u^ 5 0.005. 

The resulting steady-state solutions were checked for conservation of 
mass at each streamwise station by using the integral form of the continuity 
equation. As a percentage of the mass inflow at the upstream boundary, the 
mass loss ranged from 0.2% for the inlet and cascade flow at Re = 150 and 
Re = 0.089 to 2.3% for the inlet flow with Re = 20 and = 0.364. 
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4. Results and Discussion 

For the purpose of comparison with the solutions of Gillis and Brandt 
(1964), the u and v velocity profiles for the uniform inlet with Re = 20 
and = 0.089 are shown in figures 2 and 3, At this low Mach number, solu- 
tions to the incompressible equations of motion provide a reasonable basis for 
comparison with the inherently compressible results of the present method. 

The agreement between the two methods is excellent. The physical interpreta- 
tion of the velocity maxima observed near the inlet plane is discussed by 
Gillis and Brandt (1964). The u-velocity profile attains the well-known 
parabolic shape for an incompressible Poiseuille flow at an x/h of approxi- 
mately 2.5 (not shown) and maintains that shape to the downstream exit plane 
located at x/h = 11.0. The centerline velocity also reaches the incompressible 
asymptotic value of 1,5 u^ at x/h = 2.5. The gradual increase in the center- 
line velocity for x/h = 2,5 is due to the slight compressibility effect which 
is present even at this low Mach number. As will be shown, at higher inlet 
Mach numbers, this rate of increase in the centerline velocity is strongly 
effected and the classical parabolic shape is no longer obtained. 

Figure 4 compares the excess pressure drop at the inlet plane with the 
results of Gillis and Brandt (1964). The excess pressure drop is the pressure 
drop beyond that which would be present if the flow were a fully developed 
Poiseuille flow starting at the inlet plane, and is caused by the viscous 
interaction in the vicinity of the inlet plane. Mathematically, the excess 
pressure drop at any point (x*, y) is defined by 

aJpCx*. y) - 

x^L J 
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In the present study, 3p(x, y)/8x was evaluated numerically by using the 
specified pressure at the exit plane and the pressure at an upstream point. 

As with the velocity distribution, the agreement between the two methods is 
excellent. 

The centerline velocity distributions for the case of the uniform inlet 
= 0.089 and Re = 150 is compared to the incompressible numerical solu- 
tion of Gillis and Brandt (1964) and McDonald, Denny, and Mills (1972), and 
the analytical results of Van Dyke (1970) in figure 5. Figure 6 makes similar 
comparisons with the analytical solutions obtained by Van Dyke (1970) and the 
numerical results of Wang and Longwell (1964) (numerical solution of incom- 
pressible Navier-Stokes equations) for a cascade with Re = 150 and = 0.089 
There is good agreement between the various numerical solutions, and the 
slight deviation of the analytical results for the uniform inlet flow is 
discussed by Van Dyke (1970) . 

Velocity and static pressure distribution on the centerline and at the 
wall of a cascade with a Reynolds number of 20 and a Mach number of 0.36 is 
presented in figure 7. The behavior of the centerline velocity in this case 
is dramatically different from that which is observed for the incompressible 
case. The asymptotic value of 1.5 u^ is no longer approached; rather, the 
strong favorable pressure gradient that is required to maintain the flow also 
forces a reduction in the gas density. Conservation of mass, in turn, 
requires that the flow continuously accelerate. As a result, it is not possi- 
ble to approach the asymptotic Poiseuille type of flow when compressibility 
effects are present. The flow approaching the leading edge of the plate 
decelerates monotonically with the pressure approaching, but not quite reach- 
ing, the isentropic stagnation pressure at the leading edge. It is interesting 
to note that the pressure on the centerline in the vicinity of the inlet plane 
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is lower on the channel centerline than at y/h = 0, but then drops below 
the wall pressure inside the channel before proceeding toward a condition of 
uniform pressure, which is reached at an x/h of approximately 2. This 
behavior is similar in character to the inviscid flow around the forward por- 
tion of an airfoil, where the stagnation line static pressure first increases 
as the flow approaches the stagnation point and then drops below the free- 
stream pressure as the flow accelerates away from the leading edge. In the 
present case, the behavior is attributed to the very rapid rate of growth of 
the boundary layer near the leading edge (Gillis and Brandt 1964), which 
causes the flow near the plane of the plate to decelerate as the leading edge 
is approached. This is followed by a gradual reduction in the u-component 
of velocity as the flow turns to become parallel to the plate, which requires 
the negative normal pressure gradient observed downstream of the leading edge. 

The centerline and wall temperature distributions for this case are 
shown in figure 8. Although the centerline temperature distribution falls 
monotonically due to the expansion process, the temperature at y/h = 0.0 
first increases, recovering to approximately 60% of the isentropic stagnation 
temperature, before proceeding to fall as the flow moves down the channel. The 
difference between the adiabatic wall temperature and the centerline tempera- 
ture increases continuously, as might be expected, due to the increasing local 
Mach number. 

The effect of the location of the exit plane is also shown in figure 8, 
where temperature distributions are shown with the downstream boundary at 
x/h = 2.67 and 3.25. Where no discemable difference exists between the two 
cases, the - 2.67 results are not plotted. Only the data in the immedi- 

ate vicinity of the downstream boundary is effected by the boundary location. 
The velocity and pressure distributions for these two cases were, for the 
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scale o£ figure 7, essentially identical, including the data in the vicinity 
of the exit plane. 

Also, for the cascade with Re = 20 and M = 0.364, the u/u and v/u 

profiles are presented in figures 9 and 10, respectively. The u-velocity 

profile at the exit plane (x /h = 3.25) differs significantly from the 

e 

classical parabolic profile, in contrast with the incompressible case where 
the velocity profile is within 1% of being parabolic at x/h = 2,25 (Gillis 
and Brandt 1964). The local centerline Mach number at the exit plane is 0.94. 
Attempts to move this boundary further downstream result in numerical 
instabilities of the type described below. 

Figure 11 summarizes the results of solutions obtained for a uniform 
inlet flow at Re = 20 and several Mach numbers. The pronounced effect of 
increasing Mach number on the flow acceleration is apparent. This compressi- 
bility effect is, in fact, present even at the lowest Mach number of 0.089 
where the asymptotic incompressible value for the centerline velocity is sur- 
passed at an x/h of 2.5. At the higher inlet Mach numbers, sonic velocities 
are approached within a few channel widths. It is necessary that the down- 
stream boundary be kept upstream of the point where the flow attempts to go 
supersonic. If the downstream boundary is moved beyond this sonic point, 
numerical instability is encountered. This behavior is consistent with the 
fact that subsonic, viscous, adiabatic flow through a constant area duct can 
attain, at most, sonic velocity (see for example Shapiro 1953). To increase 
the duct length beyond this sonic point results in a "choked” flow which moves 
the sonic point downstream to the duct exit. To achieve this, a reduction in 
the mass flow rate is necessary, which is not permitted by the present 
formulation of the boundary conditions at the inlet plane. 
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5 . Summary 

The compressible viscous flow in the inlet region of a straight channel 
and the flow into a cascade of semi-infinite flat plates has been studied at 
Reynolds number of 20 and 150 and at Mach numbers ranging from 0.089 to 0.364. 
The numerical method of MacCormack (1970) is used to solve the governing 
Navier-Stokes equations utilizing a nonuniform rectangular grid. Sufficient 
boundary conditions for the solution of the governing equations are prescribed 
by including the effect of near field subsonic permeable boundaries. These 
boundary conditions are well -posed in the sense that the pressure on the 
upstream in-flow boundary is advanced in time such that accurate steady-state 
solutions are obtained. No attempt is made to relate the time-dependent 
part of the solutions to any particular physical unsteady flow. To do this for 
subsonic flow requires knowledge of conditions outside the permeable boundary. 

Comparison of the low Mach number uniform inlet and cascade results 
with the incompressible solutions of previous authors shows excellent agree- 
ment and verifies the consistency of the time-dependent boundary conditions on 
the permeable boundaries and the utility of the nonuniform grid. Results of 
computations at the higher Mach numbers for both the uniform inlet and the 
cascade are presented. The extreme sensitivity of the inlet region flow to 
increasing inlet Mach number is demonstrated by a rapid acceleration of the 
flow toward the sonic velocity. Solutions are presented in sufficient detail 
to provide basic test cases for future analytical or numerical methods which 
include the effect of compressibility on viscous flows in inlet regions. 
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Figure Captions 

FIGURE 1. Coordinate system, numerical grid, and boundary conditions. 

FIGURE 2. Streamwise velocity profiles for uniform inlet at Re = 20, = 0.089, 

, full numerical (present method) ;0> full numerical (Gillis and Brandt 

1964); X, parabolic profile. 

FIGURE 3. Transverse velocity profiles for uniform inlet at Re = 20, = 0.089. 

present method; O, Gillis and Brandt (1964). 

FIGURE 4. Excess pressure drop at inlet plane for uniform inlet at Re = 20 

and = 0.084. O, present method; A, Gillis and Brandt (1964). 

FIGURE 5. Centerline velocity for uniform inlet at Re = 150. Cj, present 

method (M^ = 0.089); , Van Dyke (1970); O, Gillis and Brandt (1964); 

, McDonald, Denny, and Mills (1972). 

FIGURE 6. Centerline velocity for cascade at Re = 150. O, present method 

^go “ 0.089; , Wang and Longwell (1964) M^ = 0; Van Dyke 

(1970) M^ = 0. 

FIGURE 7. Cascade pressure and streamwise velocity at Re = 20 and M^ = 0.364. 
D, centerline velocity; A, stagnation line velocity; 0> centerline static 
pressure; O , stagnation line and wall static pressure; , isentropic stagna- 
tion pressure. 

FIGURE 8. Cascade centerline and wall static temperature distributions at 
Re - 20 and M^ = 0.364. O , adiabatic wall temperature (x^/h = 3.25); C , 
adiabatic wall temperature (x^/h = 2.67); ED, centerline temperature (x^/h = 

3.25); A, centerline temperature (Xg/h = 2.67). 

FIGURE 9. Streamwise velocity profiles for cascade at Re = 20 and M^ = 0.364. 
— , present method; , parabolic profile. 

FIGURE 10. Transverse velocity profiles for cascade at Re = 20 and = 0.364. 
FIGURE 11. Centerline velocity distributions for uniform inlet with Re = 20. 

M = 0.089; , M = 0.196; — , M = 0.288; 

M = 0.579. 
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